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Mode  Mixity  Determinations  for  Interfacial  Cracking  in  Incompressible 
Materials  under  Plane  Strain  Conditions 

T.  C.  Miller 

Sparta,  Inc.,  4  Draco  Drive,  Edwards  Air  Force  Base,  California  93524 

Introduction 

A  frequent  site  for  the  development  of  cracks  is  along  the  interface  between  two  materials,  for  example, 
the  interface  between  the  rubber  liner  and  the  propellant  grain  in  a  solid  rocket  motor.  Cracks  along  this  interface 
experience  plane  strain  conditions  and  are  bounded  by  two  incompressible  materials.  These  conditions  lead  to 
a  simplified  analysis  so  that  complete  determination  of  the  complex  stress  intensity  factor  is  straightforward  [1,2]. 
To  study  plane  strain  interfacial  cracking  in  incompressible  material  pairs,  researchers  used  the  photoelastic  stress 
freezing  method.  Transparent  materials  which  exhibit  changes  in  the  index  of  refraction  when  stressed  are  loaded 
at  a  temperature  above  a  temperature  Tc,  known  as  the  stress  freezing  temperature.  Lowering  the  temperature 
below  Tc  while  maintaining  the  load  “locks  in”  the  stress  data  through  local  variations  in  the  index  of  refraction. 
Illumination  with  laser  light  reveals  fringe  patterns;  these  fringes  are  proportional  to  the  maximum  in-plane  shear 
stress.  Subsequent  slicing  of  the  specimen  does  not  affect  the  fringe  patterns  and  allows  for  analysis  of  the 
stresses  in  three  dimensions  [1]. 

To  analyze  plane  strain  incompressible  interfacial  fracture,  specimens  were  constructed  from  two 
materials.  The  first  was  a  plastic,  araldite,  and  the  second  was  the  same  plastic  with  the  addition  of  aluminum 
particles.  Young’s  moduli  for  these  two  materials  were  1 8.6  and  36.9  MPa,  respectively,  for  the  araldite  and  the 
araldite-aluminum,  so  that  the  two  moduli  differed  by  a  factor  of  about  two  (in  solid  rocket  motors,  the  mismatch 
of  the  two  moduli  is  between  two  and  three).  Different  mode  mixities  were  introduced  by  varying  the  crack 
orientation  in  different  specimens  while  keeping  the  loading  direction  constant.  Crack  orientations  of  T  =  0°,  15°, 
30°,  and  45°  were  considered,  where  T  =  0°  corresponds  to  mode  I  loading  (see  Fig.  1).  Analysis  of  the 
photoelastic  fringe  pattern  provided  results  for  the  magnitude  and  phase  angle  of  the  complex  stress  intensity 
factor  for  a  range  of  mode  mixities  [1]. 

Figure  1  shows  one  of  the  finite  element  meshes  used  to  model  these  specimens.  The  mesh  used  eight 
noded  plane  strain  quadrilateral  elements.  Elements  around  the  crack  tip  are  quadrilaterals  degenerated  into 
triangular  elements  with  quarter  point  nodes.  All  of  the  singularity  elements  have  the  same  node  at  the  crack  tip. 
The  notch  geometry  shown  matches  that  used  in  the  photoelastic  experiment.  Boundary  conditions  are  applied 
using  point  loads  at  the  top  of  the  aluminum  grips  (to  which  the  specimens  were  glued)  which  were  allowed  to 
rotate  freely  (as  were  the  actual  grips  used  in  the  experiment).  In  the  experiment,  an  adhesive  was  used  to  bond 
the  two  materials  together.  The  modulus  of  the  adhesive  was  similar  to  that  of  the  araldite-aluminum  (they 
differed  by  12%),  and  the  adhesive  layer  was  thin  (less  than  0.5  mm).  Because  the  layer  was  thin  and  its 
properties  were  similar  to  one  of  the  other  materials,  the  adhesive  layer  did  not  affect  the  fracture  parameters,  and 
was  not  incorporated  into  the  finite  element  models. 
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Figure  1  -  Typical  Finite  Element  Mesh  Used 
in  Numerical  Computations  (crack 
orientation  =  15  degrees). 


Mixed  Formulations  and  Incompressiblity 

Plane  strain  conditions  in  incompressible  materials  require  special  solution  techniques;  in  particular,  a 
mixed  mode  formulation  must  be  used.  The  constitutive  model  for  an  isotropic  linear  elastic  material  can  be 
expressed  as: 


Here  ay  is  the  stress  tensor,  e,j  is  the  strain  tensor,  e,/  is  the  deviatoric  strain  tensor,  and  ev  is  the 
volumetric  strain  component.  The  bulk  modulus  and  shear  modulus,  K  and  u,  characterize  the  resistance  of  the 
material  to  compression  (or  expansion)  and  shear  deformations.  However,  as  Poisson’s  ratio,  v,  approaches  14, 
the  bulk  modulus  becomes  infinite  and  the  volumetric  strain  approaches  zero.  This  makes  the  constitutive  form 
of  eq.  (1)  indeterminate,  so  that  the  stresses  require  an  additional  variable  (such  as  pressure  or  hydrostatic  stress) 


for  their  characterization  [3-5]: 


o9  =  -pby  +  2|ae/.,  p  =  -Kev  .  (2) 

When  attempting  to  model  an  incompressible  material  using  the  finite  element  method,  as  v  ->  14,  the 
diagonal  terms  of  the  local  stiffness  matrices  approach  zero,  resulting  in  an  ill-conditioned  global  stiffness  matrix. 
This  ill-conditioning  leads  to  poor  predictions  of  stress  if  a  conventional  (i.e.,  displacement  based)  finite  element 
formulation  is  used.  The  remedy  is  to  use  a  mixed  formulation,  where  the  pressure  is  a  solution  variable  as  well 
as  the  displacement  components.  The  problem  is  then  solved  by  minimizing  the  total  potential  energy,  which 
now  incorporates  the  additional  constraint  p  =  -Kev  as  a  Lagrange  multiplier  term.  Minimizing  this  revised 
potential  energy  expression  is  realized  in  the  finite  element  analysis  as  a  mixed  mode  solution  incorporating 
hybrid  elements  [3],  The  solution  variables  for  each  element  are  the  two  components  of  displacement  at  each 
node  and  three  pressure  terms  (which  allow  for  linear  pressure  variations  within  each  element)  [4,5]. 

Determination  of  Fracture  Parameters 

The  general  solution  for  a  crack  along  an  interface  between  two  isotropic  linear  elastic  materials  consists 
of  terms  related  to  an  infinite  series  of  complex  eigenvalues  [6].  In  a  region  near  the  crack  tip,  the  first  eigenvalue 
term  dominates,  so  that  without  mode  III  loading  the  asymptotic  field  equations  are  given  by  [7]: 

Op,  =  — 2^(6)  *  Im(Kr->) 2"(0))  .  (3) 

\JTkt 


Here  2^(0)  and  Epq!I(0)  are  dimensionless  mode  I  and  rpode  II  angular  functions  and  K  is  the  complex 
stress  intensity  factor,  which  is  commonly  characterized  either  by  its  real  and  imaginary  parts  or  in  a  polar  form 
(K  =  K,  +  iK2  =  Ke'V  The  coordinate  system  is  at  the  crack  tip,  with  the  x  axis  parallel  to  the  interface.  The 
two  modes  of  loading  are  coupled  in  that  their  proportions  vary  with  distance  from  the  crack  tip  for  any  given  set 
of  applied  loads.  The  bimaterial  parameter  e  is  really  the  complex  part  of  the  asymptotic  eigenvalue 
X  =  !4  +  ie,  and  is  related  to  the  second  Dundurs  parameter,  P  [2]: 


±  p  .  ^  ;_j) 

2tc  1  +  p  P,(k2  +  1)  +  P2(k,  +  1) 


(4) 


Here  p  is  the  shear  modulus,  k  =  3  -  4v  (in  plane  strain),  and  the  subscripts  index  the  two  materials.  The 
presence  of  the  complex  part  of  the  eigenvalue  causes  the  inherent  mode  coupling  and  leads  to  other 
complications  such  as  predicting  crack  face  interpenetration.  However,  substituting  v  =  14  and  k  —  3  -  4v  in  the 
equations  above  gives  e  =  P  =  0,  so  that  for  incompressible  materials  and  plane  strain  conditions,  these 
complications  vanish.  Similar  expressions  can  be  written  that  use  «£  the  complex  variables  representation.  In 
particular,  along  the  bonded  portion  of  the  interface  [8]: 


'Throughout  the  text,  K  denotes  the  complex  stress  intensity  factor,  and  K  denotes  its 
magnitude. 
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(%  +  ioxy)e=  0 


/2t vr 


Kx  +  iK2 
\j2izr 


(5) 


This  equation  shows  that  tan''[K2/K,]  and  tan"1  [(oxy/qy),=0  ]  are  equal  in  the  asymptotic  region.  The  finite  element 
data  near  the  crack  tip  can  be  used  to  determine  tan'l[K2/K]]  easily,  or  the  phase  angle  of  K,  by  extrapolating  to 
r  =  0.  To  implement  this  algorithm,  stress  components  oxy  and  are  averaged  at  the  nodes  along  the  bond  line, 
then  their  arctangent  is  plotted  versus  r  (normalized  by  the  crack  length,  a)  in  a  certain  region  near  the  crack  tip. 
In  this  work,  the  first  four  points  near  the  crack  tip  were  excluded/due  to  difficulties  with  adequately  modeling 
the  high  gradients  very  near  4$  the  crack  tip.  Also,  points  with  r/a  >  1  were  excluded;  far  from  the  crack  the 
stresses  approximate  those  of  a  plate  subjected  to  combined  tension  and  shear  loadingjor  approach  a  traction  free 
surface  condition  (as  the  free  surface  is  approached).  The  remaining  data  is  fit  with  a  cubic  polynomial,  and  the 
limit  is  taken  as  r  -  0,  so  that  the  constant  term  in  the  polynomial  is  an  estimate  of1?: 


Y  =  I/w(r- 0)  [(%/%)0=o]  a  Lim(r/a- 0)  i|r(r/a)  .  (6) 

Here  Y  is  the  complex  stress  intensity  factor  phase  angle  and  i|i  is  the  fitting  function  for 
tan'1  [(o^Oyy)  e=o]-  111  practice  both  the  fitting  function  and  its  domain  can  be  arbitrarily  chosen  provided  that  they 
adequately  represent  the  data  near  the  crack  tip. 

Magnitudes  can  also  be  determined  easily  when  e  =  0.  The  magnitude  of  K  is  related  to  the  energy  release 
rate  [8]: 


J  =  G  = 


=  =  y[=  +  =] 
E'  2  Ex  E2 


(7) 


The  J  integral  can  be  easily  calculated  using  techniques  such  as  virtual  crack  extension  or  the  domain 
integral  method,  mi  K  can  then  be  determined.  Finite  element  data  can  be  used  to  characterize  the  complex 
stress  intensity  factor  completely  by  (a)  simple  curve  fitting  and  extrapolation  of  bond  line  traction  data  to 
determine  the  phase  angle  and  (b)  deriving  the  magnitude  from  J  values  and  the  effective  plane  strain  modulus, 
E*. 


Results 

Figure  2  shows  oyy  contour  plots  for  each  of  the  four  geometries  considered  (in  these  plots,  each  specimen 
has  been  subjected  to  a  remotely  applied  tensile  stress  of  9.65  kPa).  The  mode  mixity  for  the  crack  increases 
approximately  linearly  with  the  crack  orientation  angle  T,  so  that  the  contour  plots  show  the  variance  of  oyy  with 
mode  mixity.  A  comparison  of  the  four  contour  plots  shows  that  the  size  of  the  contours  decreases  as  the 
loading  changes  from  mode  I  to  mode  II.  Figure  3  shows  the  commensurate  increase  in  the  size  of  the  oxy 
contours.  The  invariance  in  the  shapes  of  the  contours  shown  in  these  two  figures  suggests  that  the  asymmetries 
shown  here  are  caused  by  the  mismatch  of  the  two  elastic  moduli  rather  than  changes  in  the  mode  mixity.  This 
shows  that  experimental  methods  that  use  fringes  proportional  to  either  o^  or  oxy  cannot  be  used  to  determine 
the  mode  mixity.  The  shape  of  the  maximum  in-plane  shear  stress  contours  depends  strongly  on  the  mode  mixity, 


so  that  experimental  methods  that  exploit  this  (such  as  the  photoelastic  stress  freezing  method)  can  be  used  for 
mode  mixity  determination. 

For  the  contour  plots  shown,  the  applied  stresses  were  all  identical,  but  in  the  photoelastic  experiment  the 
loads  wore  varied.  Experimental  results  presented  here  are  for  slices  taken  from  the  mid-thickness  of  the 
specimens,  so  that  plane  strain  conditions  were  being  measured  (the  finite  element  models  use  plane  strain 
elements  to  model  these  slices).  Table  1  gives  a  comparison  of  the  experimental  and  computational  results. 
Residual  and  machining  stresses  induced  during  specimen  fabrication  caused  the  high  experimental  stress 
intensity  magnitudes  shown  for  the  T  =  0°  and  15°  specimens.  These  stresses  cause  "no  load"  fringes  that  lie 
parallel  to  the  bond  line.  When  T  =  0°,  these  residual  stress  fringes  affect  the  gradient  of  the  fringes  in  the  y 
direction,  which  is  used  to  calculate  K,  but  this  effect  diminishes  as  the  crack  orientation  angle  increases.  The 
finite  element  model  assumes  an  initial  state  free  from  these  residual  stresses,  so  that  the  computational 
calculations  for  K  are  perhaps  more  accurate  than  the  experimental  results  for  the  lower  values  of  T.  However, 
additional  experimental  testing  can  be  used  to  adjust  for  the  effects  of  residual  and  machining  stresses  if 
necessary.  Table  1  also  shows  the  phase  angle  comparisons.  The  computational  and  experimental  results  agree 
within  5°  for  each  specimen  tested.  This  suggests  that  the  finite  element  models  can  be  used  to  characterize  K 
for  other  similar  geometries,  such  as  those  found  in  solid  rocket  motors. 


crack  orientation 
[deg] 

nominal 

stress 

fkPa] 

magnitudes  [kPa  m1/2] 

computational  experimental 

phase  angles  [deg.] 

computational  experimental 

0 

60.3 

15.8 

19.0 

4.0 

0.0 

15 

96.5 

24.0 

30.2 

11.5 

7.8 

30 

96.5 

20.8 

20.3 

19.5 

17.1 

45 

48.3 

7.2 

8.3 

25.8 

30.0 

Conclusions 


When  a  crack  experiences  plane  strain  conditions  and  lies  along  the  interface  between  incompressible 
materials,  the  near  tip  field  equations  are  simplified  by  the  vanishing  of  the  bimaterial  parameter,  e.  Finite 
element  models  that  use  a  mixed  formulation  can  be  used  to  characterize  the  complex  stress  intensity  factor  of 
these  cracks  completely.  This  is  done  by  using  J  integral  calculations  and  a  regression  of  the  ratio  of  bond  line 
traction  components  to  r  =  0. 
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